{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train a Model to Predict Formation Energy using the OQMD\n",
    "This notebook recreates a 2016 paper by [Ward et al.](https://www.nature.com/articles/npjcompumats201628) on predicting the formation enthalpy of materials based on their composition. We will use the [Materials Data Facility](http://materialsdatafacility.org) to retrieve a training set from the the [OQMD](http://oqmd.org), compute features based on the composition of each entry, and then train a random forest model.\n",
    "\n",
    "This example was last updated on 11/15/18 for Matminer v.0.4.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matminer.data_retrieval import retrieve_MDF\n",
    "from matminer.featurizers.base import MultipleFeaturizer\n",
    "from matminer.featurizers import composition as cf\n",
    "from matminer.featurizers.conversions import StrToComposition\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pickle as pkl\n",
    "from sklearn import metrics\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.model_selection import cross_val_score, cross_val_predict, GridSearchCV, ShuffleSplit, KFold"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Settings to change"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "quick_demo = True # Whether to run an faster version of this demo. \n",
    "# The full OQMD model takes about a hour to test and ~8GB of RAM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Training Set\n",
    "Ward _et al._ trained their machine learning models on the formation enthalpies of crystalline compounds form the [OQMD](oqmd.org). Here, we extract the data using the copy of the OQMD available through the MDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download the Data\n",
    "We first create a `Forge` instance, which simplifies performing search queries against the MDF."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is to create a tool for reading from the MDF's search index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "mdf = retrieve_MDF.MDFDataRetrieval(anonymous=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we assemble a query that gets only the converged static calculations from the OQMD. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "query_string = 'mdf.source_name:oqmd AND (oqmd.configuration:static OR '\\\n",
    "    'oqmd.configuration:standard) AND dft.converged:True'\n",
    "if quick_demo:\n",
    "    query_string += \" AND mdf.scroll_id:<10000\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = mdf.get_data(query_string, unwind_arrays=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tool creates a DataFrame object with the metadata for each entry in the OQMD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>crystal_structure.cross_reference.icsd</th>\n",
       "      <th>crystal_structure.number_of_atoms</th>\n",
       "      <th>crystal_structure.space_group_number</th>\n",
       "      <th>crystal_structure.volume</th>\n",
       "      <th>dft.converged</th>\n",
       "      <th>dft.cutoff_energy</th>\n",
       "      <th>dft.exchange_correlation_functional</th>\n",
       "      <th>files</th>\n",
       "      <th>material.composition</th>\n",
       "      <th>material.elements</th>\n",
       "      <th>...</th>\n",
       "      <th>oqmd.delta_e.units</th>\n",
       "      <th>oqmd.delta_e.value</th>\n",
       "      <th>oqmd.magnetic_moment.units</th>\n",
       "      <th>oqmd.magnetic_moment.value</th>\n",
       "      <th>oqmd.stability.units</th>\n",
       "      <th>oqmd.stability.value</th>\n",
       "      <th>oqmd.total_energy.units</th>\n",
       "      <th>oqmd.total_energy.value</th>\n",
       "      <th>oqmd.volume_pa.units</th>\n",
       "      <th>oqmd.volume_pa.value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>165739.0</td>\n",
       "      <td>4</td>\n",
       "      <td>225</td>\n",
       "      <td>58.7515</td>\n",
       "      <td>True</td>\n",
       "      <td>520.0</td>\n",
       "      <td>PBE</td>\n",
       "      <td>[{'data_type': 'ASCII text, with very long lin...</td>\n",
       "      <td>Al1Cu1S2</td>\n",
       "      <td>[Al, Cu, S]</td>\n",
       "      <td>...</td>\n",
       "      <td>eV/atom</td>\n",
       "      <td>-0.066756</td>\n",
       "      <td>bohr/atom</td>\n",
       "      <td>NaN</td>\n",
       "      <td>eV/atom</td>\n",
       "      <td>0.893649</td>\n",
       "      <td>eV/atom</td>\n",
       "      <td>-3.851622</td>\n",
       "      <td>angstrom^3/atom</td>\n",
       "      <td>14.6879</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>75725.0</td>\n",
       "      <td>12</td>\n",
       "      <td>123</td>\n",
       "      <td>180.0000</td>\n",
       "      <td>True</td>\n",
       "      <td>520.0</td>\n",
       "      <td>PBE</td>\n",
       "      <td>[{'data_type': 'ASCII text, with very long lin...</td>\n",
       "      <td>Ba2Ca1Cu2Hg1O6</td>\n",
       "      <td>[Ba, Ca, Cu, Hg, O]</td>\n",
       "      <td>...</td>\n",
       "      <td>eV/atom</td>\n",
       "      <td>-1.851396</td>\n",
       "      <td>bohr/atom</td>\n",
       "      <td>NaN</td>\n",
       "      <td>eV/atom</td>\n",
       "      <td>0.035518</td>\n",
       "      <td>eV/atom</td>\n",
       "      <td>-5.004761</td>\n",
       "      <td>angstrom^3/atom</td>\n",
       "      <td>15.0000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2 rows × 31 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   crystal_structure.cross_reference.icsd  crystal_structure.number_of_atoms  \\\n",
       "0                                165739.0                                  4   \n",
       "1                                 75725.0                                 12   \n",
       "\n",
       "   crystal_structure.space_group_number  crystal_structure.volume  \\\n",
       "0                                   225                   58.7515   \n",
       "1                                   123                  180.0000   \n",
       "\n",
       "   dft.converged  dft.cutoff_energy dft.exchange_correlation_functional  \\\n",
       "0           True              520.0                                 PBE   \n",
       "1           True              520.0                                 PBE   \n",
       "\n",
       "                                               files material.composition  \\\n",
       "0  [{'data_type': 'ASCII text, with very long lin...             Al1Cu1S2   \n",
       "1  [{'data_type': 'ASCII text, with very long lin...       Ba2Ca1Cu2Hg1O6   \n",
       "\n",
       "     material.elements          ...          oqmd.delta_e.units  \\\n",
       "0          [Al, Cu, S]          ...                     eV/atom   \n",
       "1  [Ba, Ca, Cu, Hg, O]          ...                     eV/atom   \n",
       "\n",
       "  oqmd.delta_e.value oqmd.magnetic_moment.units oqmd.magnetic_moment.value  \\\n",
       "0          -0.066756                  bohr/atom                        NaN   \n",
       "1          -1.851396                  bohr/atom                        NaN   \n",
       "\n",
       "   oqmd.stability.units oqmd.stability.value oqmd.total_energy.units  \\\n",
       "0               eV/atom             0.893649                 eV/atom   \n",
       "1               eV/atom             0.035518                 eV/atom   \n",
       "\n",
       "   oqmd.total_energy.value oqmd.volume_pa.units  oqmd.volume_pa.value  \n",
       "0                -3.851622      angstrom^3/atom               14.6879  \n",
       "1                -5.004761      angstrom^3/atom               15.0000  \n",
       "\n",
       "[2 rows x 31 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We only need two columns: `delta_e` and `material.composition`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data[['oqmd.delta_e.value', 'material.composition']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Renaming the columns to make the rest of the code more succinct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data.rename(columns={'oqmd.delta_e.value': 'delta_e', 'material.composition':'composition'})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compile the Training Set\n",
    "Our next step is to get only the lowest-energy entry for each composition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4561f6c9af914aefb501a196b744f971",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='StrToComposition', max=4849), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "data = StrToComposition(target_col_id='composition_obj').featurize_dataframe(data, 'composition')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create shortcuts for our input and output columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove compounds w/o a `delta_e` measurement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k in ['delta_e']:\n",
    "    data[k] = pd.to_numeric(data[k])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Removed 409/4849 entries\n"
     ]
    }
   ],
   "source": [
    "original_count = len(data)\n",
    "data = data[~ data['delta_e'].isnull()]\n",
    "print('Removed %d/%d entries'%(original_count - len(data), original_count))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get only the groundstate and each composition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Removed 24/4440 entries\n",
      "CPU times: user 643 ms, sys: 8.83 ms, total: 652 ms\n",
      "Wall time: 742 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "original_count = len(data)\n",
    "data['composition'] = data['composition_obj'].apply(lambda x: x.reduced_formula)\n",
    "data.sort_values('delta_e', ascending=True, inplace=True)\n",
    "data.drop_duplicates('composition', keep='first', inplace=True)\n",
    "print('Removed %d/%d entries'%(original_count - len(data), original_count))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Removed 1/4416 entries\n"
     ]
    }
   ],
   "source": [
    "original_count = len(data)\n",
    "data = data[np.logical_and(data['delta_e'] >= -20, data['delta_e'] <= 5)]\n",
    "print('Removed %d/%d entries'%(original_count - len(data), original_count))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build an ML model\n",
    "In this part of the notebook, we build a ML model using [scikit-learn](http://scikit-learn.org/stable/) and evaluate its performance using cross-validation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 1: Compute Representation\n",
    "The first step in building a ML model is to convert the raw materials data (here: the composition) into the required input for an ML model: a finite list of quantitative attributes. In this example, we use the \"general-purpose\" attributes of [Ward *et al* 2016](https://www.nature.com/articles/npjcompumats201628)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_calculators = MultipleFeaturizer([cf.Stoichiometry(), cf.ElementProperty.from_preset(\"magpie\"),\n",
    "                                          cf.ValenceOrbital(props=['avg']), cf.IonProperty(fast=True)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the feature names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_labels = feature_calculators.feature_labels()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fa9181e6b9c442a39cbc6bd7a9a3893a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='MultipleFeaturizer', max=4415), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "CPU times: user 812 ms, sys: 176 ms, total: 988 ms\n",
      "Wall time: 20.8 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "data = feature_calculators.featurize_dataframe(data, col_id='composition_obj');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generated 145 features\n",
      "Training set size: 4415x145\n"
     ]
    }
   ],
   "source": [
    "print('Generated %d features'%len(feature_labels))\n",
    "print('Training set size:', 'x'.join([str(x) for x in data[feature_labels].shape]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove entries with `NaN` or `infinite` features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Removed 0/4415 entries\n"
     ]
    }
   ],
   "source": [
    "original_count = len(data)\n",
    "data = data[~ data[feature_labels].isnull().any(axis=1)]\n",
    "print('Removed %d/%d entries'%(original_count - len(data), original_count))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 2: Tuning Model Hyper-parameters\n",
    "For brevity, we will only consider one ML algorithm in this example: [random forest](http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html#sklearn.ensemble.RandomForestRegressor). The \"random forest\" algorithm works by training many different decision tree models, where each is trained on a different subset of the dataset . Here, we tune one of the major parameters of the algorithm: the number features considered at each split in each decision tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = GridSearchCV(RandomForestRegressor(n_estimators=20 if quick_demo else 150, n_jobs=-1),\n",
    "                     param_grid=dict(max_features=range(8,15)),\n",
    "                     scoring='neg_mean_squared_error',cv=ShuffleSplit(n_splits=1, test_size=0.1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GridSearchCV(cv=ShuffleSplit(n_splits=1, random_state=None, test_size=0.1, train_size=None),\n",
       "       error_score='raise-deprecating',\n",
       "       estimator=RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,\n",
       "           max_features='auto', max_leaf_nodes=None,\n",
       "           min_impurity_decrease=0.0, min_impurity_split=None,\n",
       "           min_samples_leaf=1, min_samples_split=2,\n",
       "           min_weight_fraction_leaf=0.0, n_estimators=20, n_jobs=-1,\n",
       "           oob_score=False, random_state=None, verbose=0, warm_start=False),\n",
       "       fit_params=None, iid='warn', n_jobs=None,\n",
       "       param_grid={'max_features': range(8, 15)}, pre_dispatch='2*n_jobs',\n",
       "       refit=True, return_train_score='warn',\n",
       "       scoring='neg_mean_squared_error', verbose=0)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.fit(data[feature_labels], data['delta_e'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the tuning results. This shows the CV score as a function of the parameter we tuned \"max features\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.08189797053382962"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.best_score_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'RMSE (eV/atom)')"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAHQdJREFUeJzt3X+YHnV97vH3TRIgwimxEKRJCGCMsQEi4Aq1iD8QJIolAUUTpKJSESXallMqOXBoxUtRYvWUGivURmwVIgJiTkVDD8UfoHiyIcEQIBB+KElUosKJkQBZuM8fM4tPlt2dJ5vMPrvP3q/r2muf+c7MM5+5xNw78535fmWbiIiI/uzS6gIiImLoS1hERESlhEVERFRKWERERKWERUREVEpYREREpYRFRERUqjUsJM2UtEbSWknn97L+bEmrJK2UdKuk6T3WT5a0WdLf1FlnRET0T3W9lCdpFHAfcDywDlgGzLV9d8M2f2B7U/n5JOCDtmc2rL8WMPBj25+updCIiKg0usbvPhJYa/tBAEmLgVnAc2HRHRSlPSiCgXL72cBDwO+aOdg+++zjAw88cMerjogYQZYvX/4r2+OrtqszLCYCjzQsrwOO6rmRpHOAc4FdgWPLtj2Bj1BclfR5C0rSWcBZAJMnT6azs3Nn1R4RMSJI+mkz27W8g9v2QttTKMLhwrL574HP2t5cse8Vtjtsd4wfXxmMERExQHVeWawH9m9YnlS29WUx8M/l56OAt0m6FBgHPCvpSdufq6XSiIjoV51hsQyYKukgipCYA5zWuIGkqbbvLxdPBO4HsH1MwzZ/D2xOUEREtE5tYWG7S9I8YCkwClhke7Wki4FO20uAeZKOA7YCjwFn1FVPREQMXG2Pzg62jo4Op4M7ImL7SFpuu6Nqu5Z3cEdExNCXsIiIiEoJi4iIqJSwiIiISgmLiIiolLCIiIhKCYuIiKiUsIiIiEoJi4iIqJSwiIiISgmLiIiolLCIiIhKCYuIiKiUsIiIiEoJi4iIqJSwiIiISgmLiIiolLCIiIhKCYuIiKiUsIiIiEoJi4iIqFRrWEiaKWmNpLWSzu9l/dmSVklaKelWSdPL9iPLtpWS7pR0cp11RkRE/2oLC0mjgIXAm4DpwNzuMGhwle1DbR8GXAp8pmy/C+go22cCl0saXVetERHRvzqvLI4E1tp+0PbTwGJgVuMGtjc1LO4BuGx/wnZX2b57d3tERLRGnX+tTwQeaVheBxzVcyNJ5wDnArsCxza0HwUsAg4A/rwhPCIiYpC1vIPb9kLbU4CPABc2tP/Y9sHAK4H5knbvua+ksyR1SurcuHHj4BUdETHC1BkW64H9G5YnlW19WQzM7tlo+x5gM3BIL+uusN1hu2P8+PE7WG5ERPSlzrBYBkyVdJCkXYE5wJLGDSRNbVg8Ebi/bD+ou0Nb0gHAy4CHa6w1IiL6UVufhe0uSfOApcAoYJHt1ZIuBjptLwHmSToO2Ao8BpxR7v5q4HxJW4FngQ/a/lVdtUZERP9kt8eDRh0dHe7s7Gx1GRERw4qk5bY7qrZreQd3REQMfQmLiIiolLCIiIhKCYuIiKiUsIiIiEoJi4iIqJSwiIiISgmLiIiolLCIiIhKCYuIiKiUsIiIiEoJi4iIqJR5rSMihqkbVqxnwdI1bHh8CxPGjeW8E6Yx+/CJtRwrYRERMQzdsGI9869fxZatzwCw/vEtzL9+FUAtgZHbUBERw9CCpWueC4puW7Y+w4Kla2o5XsIiImIY2vD4lu1q31EJi4iIYWjCuLHb1b6jEhYRMaLcsGI9R3/yvzjo/G9x9Cf/ixtWrG91SQNy3gnTGDtm1DZtY8eM4rwTptVyvHRwR8SIMdidwnXqrjdPQ0VE7GT9dQoPt7CAIjAGq+7choqIEWOwO4XbScIiIkaMwe4Ubie1hoWkmZLWSFor6fxe1p8taZWklZJulTS9bD9e0vJy3XJJx9ZZZ0SMDIPdKdxOauuzkDQKWAgcD6wDlklaYvvuhs2usv2FcvuTgM8AM4FfAX9me4OkQ4ClwPC7oRgRQ8pgdwq3kzo7uI8E1tp+EEDSYmAW8FxY2N7UsP0egMv2FQ3tq4Gxknaz/VSN9UbECDCYncLtpM6wmAg80rC8Djiq50aSzgHOBXYFervd9Fbgjt6CQtJZwFkAkydP3gklR0REb1rewW17oe0pwEeACxvXSToY+BTw/j72vcJ2h+2O8ePH119sRMQIVWdYrAf2b1ieVLb1ZTEwu3tB0iTgG8C7bD9QS4UREdGUOsNiGTBV0kGSdgXmAEsaN5A0tWHxROD+sn0c8C3gfNu31VhjREQ0obawsN0FzKN4kuke4BrbqyVdXD75BDBP0mpJKyn6Lc7obgdeAlxUPla7UtK+ddUaERH9k+1W17BTdHR0uLOzs9VlREQMK5KW2+6o2m7Ejw01mNMSRkQMVyM6LNppBMqIiDptV5+FpD3KN7PbwmBPSxgRMVz1GxaSdpF0mqRvSXoUuBf4uaS7JS2Q9JLBKbMeGYEyIqI5VVcWtwBTgPnAfrb3t70v8GrgduBTkk6vucbaZATKiIjmVPVZHGd7a89G278BrgOukzSmlsoGwXknTNumzwIyAmVERG/6DYvGoJD0Qoo3skc3rL+jtzAZLjICZUREc5p6GkrSx4B3Aw9Qjgxb/h7280xkBMqIiGrNPjr7dmCK7afrLCYiIoamZh+dvQsYV2chERExdDV7ZXEJsELSXcBz80rYPqnvXSIiol00GxZfpphXYhXwbH3lRETEUNRsWDxh+7JaK4mIiCGr2bD4gaRLKOajaLwNdUctVcWAZFDEiKhLs2FxePn7Txra2uLR2XaRQREjok5NhYXt19ddSOyY/gZFTFhExI5q6tFZSXtJ+oykzvLnHyTtVXdx0bwMihgRdWr2PYtFwG8pXs57O7AJ+FJdRcX2y6CIEVGnZsNiiu2/s/1g+fNR4MV1Fhbb57wTpjF2zLZTjWRQxIjYWZoNiy2SXt29IOloIPc3hpDZh0/kklMOZeK4sQiYOG4sl5xyaPorImKnaPZpqLOBf2vop3gMOKOekmKgMihiRNSl2SuLTbZfDswAZtg+nKIPo1+SZkpaI2mtpPN7WX+2pFWSVkq6VdL0sn1vSbdI2izpc9tzQhERsfM1GxbXAdjeZHtT2XZtfzuUc3UvBN4ETAfmdodBg6tsH2r7MOBS4DNl+5PA/wT+psn6IiKiRv3ehpL0MuBgYC9JpzSs+gNg94rvPhJYa/vB8rsWA7OAu7s3aAgegD0o58qw/Tvg1uE+x3dERLuo6rOYBryFYnjyP2to/y3wvop9JwKPNCyvA47quZGkc4BzgV3JG+EREUNS1bSq3wS+KelVtn9URwG2FwILJZ0GXMh2dJxLOgs4C2Dy5Ml1lBcRETT/NNSK8grgYBpuP9l+bz/7rKeYs7vbpLKtL4uBf26ynu7jXwFcAdDR0eGKzSMiYoCa7eD+d2A/4ATgexT/8Fc9DbUMmCrpIEm7AnMoRq19jqSpDYsnAvc3WU9ERAyiZq8sXmL7VEmzbH9Z0lXAD/rbwXaXpHnAUmAUsMj2akkXA522lwDzJB0HbKXHuxuSHqboSN9V0mzgjbbv7nmciIioX7NhsbX8/bikQ4BfAPtW7WT7RuDGHm0XNXz+y372PbDJ2iIiombNhsUVkl5I8e7DEmDP8nNERIwAVe9ZvAq43fYXy6bvkQEEIyJGnKoO7ncByyUtlvRuSfsNRlERETG0VL1n8QF47k3uNwFXloMJ3gJ8B7jN9jP9fEVERLSBph6dtX2v7c/anknxlvWtwKnAj+ssLiIihoaqPosbgauAG2xvBrC9heIJpxv72zciItpH1ZXF5RQvyz0o6RpJJ5cv2EVExAjSb1jY/qbtucCBFMOUvwv4maQvSTp+EOqLiIghoNk+iydsf832ycAbgcMoOrgjImIEaCosJL1I0ock3QbcQDGExxG1VhYREUNGVQf3+4C5FPNaXAecZ/uHg1FYREQMHVXDfbwKuAS42fazg1BPREQMQVUv5b0XQIXTgRfbvljSZGA/2/93MIqMiIjWanY+i89TXGXMLZd/CyyspaKIiBhymh119ijbR0haAWD7sbxvERExcjR7ZbFV0ijAAJLGA+nDiIgYIZoNi8uAbwD7Svo4xdhQn6itqoiIGFKaug1l+6uSlgNvAATMtn1PrZVFRMSQUfWexZ4NAwjeC9zb3zYREdGeqm5DfVPSP0h6jaQ9uhslvVjSmZKWAjPrLTEiIlqt6j2LN0h6M/B+4GhJfwhsBdYA3wLOsP2L+suMiIhWquyzsJ25KyIiRrhmn4YaEEkzJa2RtFbS+b2sP1vSKkkrJd0qaXrDuvnlfmsknVBnnRER0b/awqJ8L2Mhxdzd04G5jWFQusr2obYPAy4FPlPuOx2YAxxM0Sfy+fL7IiKiBeq8sjgSWGv7QdtPA4uBWY0b2N7UsLgH5Ut/5XaLbT9l+yFgbfl9ERHRAv2GhaRjGz4f1GPdKRXfPRF4pGF5XdnW8xjnSHqA4sriw9uzb0REDI6qK4tPN3y+rse6C3dGAbYX2p4CfGR7v1PSWZI6JXVu3LhxZ5QTERG9qAoL9fG5t+We1gP7NyxPKtv6shiYvT372r7CdoftjvHjx1eUExERA1UVFu7jc2/LPS0Dpko6qByhdg6wpHEDSVMbFk8E7i8/LwHmSNqtvP01FcjcGRERLVL1nsWLJS2huIro/ky5fFDfu4HtLknzKObrHgUssr1a0sVAp+0lwDxJx1G86PcYcEa572pJ1wB3A13AObafGdgpRkTEjpLd9wWCpNf2t7Pt7+30igaoo6PDnZ2drS4jImJYkbTcdkfVdlXDfWwTBpLGAIcA620/umMlRkTEcFH16OwXJB1cft4LuBP4N2CFpLn97RsREe2jqoP7GNury8/vAe6zfSjwCuBva60sIiKGjKqweLrh8/HADQAZaTYiYmSpCovHJb1F0uHA0cB3ACSNBsbWXVxERAwNVY/Ovp9i/u39gL9quKJ4A8V8FhERMQJUPQ11H73MhGd7KcX7ExERMQJUzcF9WX/rbX+4v/UREdEeqm5DnQ3cBVwDbKB6PKiIiGhDVWHxR8CpwDsoht34GnCt7cfrLiwiIoaOfp+Gsv1r21+w/XqK9yzGAXdL+vNBqS4iIoaEqisLACQdAcyleNfi28DyOouKiIihpaqD+2KKocPvoZhvYr7trsEoLCIiho6qK4sLgYeAl5c/n5AERUe3bc+ot7yIiBgKqsKi3zkrIiJiZKh6Ke+nvbVL2oWiD6PX9RE74oYV61mwdA0bHt/ChHFjOe+Eacw+fGKry4oY0aqGKP8DSfMlfU7SG1X4EPAg8PbBKTFGkhtWrGf+9atY//gWDKx/fAvzr1/FDSv6m749IupWNZDgvwPTgFXAXwC3AG8DZtueVXNtMQItWLqGLVu3nUF3y9ZnWLB0TYsqighoYg7ucv4KJH0R+Dkw2faTtVcWI9KGx7dsV3tEDI6qK4ut3R9sPwOsS1BEnSaM633k+77aI2JwVIXFyyVtKn9+C8zo/ixp02AUGCPLeSdMY+yYUdu0jR0zivNOmNaiiiICqp+GGtXf+iqSZgL/CIwCvmj7kz3Wn0vRF9IFbATe2/0ElqRPUbwQCPAx21/bkVpieOh+6qldnobKk13RLpoa7mMgJI0CFlIMEbIOWCZpie27GzZbAXTYfkLSB4BLgXdIOhE4AjgM2A34rqRv287VzAgw+/CJbfEPaveTXd0d9t1PdgFtcX4xslTdhtoRRwJrbT9o+2mK4UK2eYLK9i22nygXbwcmlZ+nA9+33WX7d8BP6GUSpoihLE92RTupMywmAo80LK8r2/pyJsUghQB3AjMlvUDSPsDrgf1rqTKiJnmyK9pJbbehtoek04EO4LUAtm+S9ErghxR9GT8Cnullv7OAswAmT548aPVGNGPCuLGs7yUY8mRXDEd1XlmsZ9urgUll2zYkHQdcAJxk+6nudtsft32Y7eMpBi68r+e+tq+w3WG7Y/z48Tv9BCJ2RJ7sinZS55XFMmCqpIMoQmIOcFrjBpIOBy4HZtp+tKF9FDDO9q8lzQBmADfVWGvETtduT3bFyFZbWNjukjQPWErx6Owi26vLOTI6bS8BFgB7Al8vhz7/me2TgDHAD8q2TcDpmUcjhqN2ebIrotY+C9s3Ajf2aLuo4fNxfez3JMUTURERMQTU2WcRERFtImERERGVEhYREVEpYRERzdm8Ge67r/gdI07CIiL619XFA3Pew5Mv3JvNhxzGky/cmwfmvAe68oDiSDIk3uCOiKHrgdPfx4Trrmb3rqefa5tw3dU8AExZ/KXWFRaDKlcWEdG3zZuZeN1VjO16apvmsV1PMfH6q3NLagRJWERE3zZsoEu9T2vTpV1gw4ZBLihaJWEREX2bMIHRft4YngCM9rMwYcIgFxStkrCIiL7tuSfr33oaW0bvtk3zltG7sf6UubDnni0qLAZbOrgjol9TvvIvPABMvP5qurQLo/0sG06Zy5Sv/EurS4tBlLCIiP6NHl089bT5n4o+igkTmJIrihEnYRERzdlzT3jpS1tdRbRI+iwiIqJSwiIiIiolLCIiolLCIiIiKiUsIiKiUsIiIiIqJSwiIqJSwiIiIiolLCIiolKtYSFppqQ1ktZKOr+X9edKulvSTyTdLOmAhnWXSlot6R5Jl0lSnbVGRETfagsLSaOAhcCbgOnAXEnTe2y2AuiwPQO4Fri03PdPgaOBGcAhwCuB19ZVa0RE9K/OK4sjgbW2H7T9NLAYmNW4ge1bbD9RLt4OTOpeBewO7ArsBowBflljrRER0Y86w2Ii8EjD8rqyrS9nAt8GsP0j4Bbg5+XPUtv31FRnRERUGBId3JJOBzqABeXyS4A/prjSmAgcK+mYXvY7S1KnpM6NGzcOZskRESNKnWGxHti/YXlS2bYNSccBFwAn2e6eFf5k4Hbbm21vprjieFXPfW1fYbvDdsf48eN3+glEREShzrBYBkyVdJCkXYE5wJLGDSQdDlxOERSPNqz6GfBaSaMljaHo3M5tqIiIFqktLGx3AfOApRT/0F9je7WkiyWdVG62ANgT+LqklZK6w+Ra4AFgFXAncKft/11XrRER0T/ZbnUNO0VHR4c7OztbXUZExLAiabntjqrthkQHd0REDG0Ji4iIqJSwiIiISgmLiIiolLCIiIhKCYuIiKiUsIiIiEoJi4iIqJSwiIiISgmLiIiolLCIiIhKCYuIiKiUsIiIiEoJi4iIqJSwiIiISgmLiIiolLCIiIhKCYuIiKiUsIiIiEoJi4iIqJSwiIiISgmLiIioVGtYSJopaY2ktZLO72X9uZLulvQTSTdLOqBsf72klQ0/T0qaXWetERHRt9rCQtIoYCHwJmA6MFfS9B6brQA6bM8ArgUuBbB9i+3DbB8GHAs8AdxUV60REdG/Oq8sjgTW2n7Q9tPAYmBW4wZlKDxRLt4OTOrle94GfLthu4iIGGR1hsVE4JGG5XVlW1/OBL7dS/sc4OredpB0lqROSZ0bN24ccKEREdG/IdHBLel0oANY0KP9j4BDgaW97Wf7CtsdtjvGjx9ff6ERESPU6Bq/ez2wf8PypLJtG5KOAy4AXmv7qR6r3w58w/bW2qqMiIhKsl3PF0ujgfuAN1CExDLgNNurG7Y5nKJje6bt+3v5jtuB+bZvaeJ4G4Gf7kDJ+wC/2oH9h5J2OZd2OQ9on3Npl/OAnEu3A2xX3pqpLSwAJL0Z+F/AKGCR7Y9LuhjotL1E0v+huM3083KXn9k+qdz3QOA2YH/bz9ZW5O9r7bTdUfdxBkO7nEu7nAe0z7m0y3lAzmV71XkbCts3Ajf2aLuo4fNx/ez7MP13iEdExCAZEh3cERExtCUsfu+KVhewE7XLubTLeUD7nEu7nAfkXLZLrX0WERHRHnJlERERlRIWgKS/lrRa0l2Srpa0e6trGghJf1mew2pJf9XqeraHpEWSHpV0V0PbH0r6T0n3l79f2Moam9XHuZxa/u/yrKRh8QROH+exQNK95eCf35A0rpU1NquPc/lYeR4rJd0kaUIra2xGb+fRsO6/S7Kkfeo49ogPC0kTgQ9TDGh4CMVjvnNaW9X2k3QI8D6KMbleDrxF0ktaW9V2uRKY2aPtfOBm21OBm8vl4eBKnn8udwGnAN8f9GoG7kqefx7/CRxSDv55HzB/sIsaoCt5/rkssD2jHLD0P4CLnrfX0HMlzz8PJO0PvBH4WV0HHvFhURoNjC1fJHwBsKHF9QzEHwM/tv2E7S7gexT/OA0Ltr8P/KZH8yzgy+XnLwPDYpj63s7F9j2217SopAHp4zxuKv/7gr4H/xxy+jiXTQ2LewBDvgO3j/+fAHwW+FtqPIcRHxa21wOfpkjknwP/z/ZwHA79LuAYSXtLegHwZrYdbmU4epHt7hc2fwG8qJXFxPO8l94H/xw2JH1c0iPAOxkeVxbPI2kWsN72nXUeZ8SHRXkffBZwEDAB2KMc2HBYsX0P8CmKeT++A6wEnmlpUTuRi8f2hvxffiOFpAuALuCrra5lR9i+wPb+FOcxr9X1bK/yD8P/wSAE3YgPC+A44CHbG8sBC68H/rTFNQ2I7X+1/QrbrwEeo7inPJz9shx5uHsE4kdbXE8Akt4NvAV4p9vn2fuvAm9tdREDMIXiD907JT1McVvwDkn77ewDJSyK209/IukFkkQx8OE9La5pQCTtW/6eTNFfcVVrK9phS4Azys9nAN9sYS1BMVUyxb3xk4b7hGSSpjYszgLubVUtA2V7le19bR9o+0CKeYOOsP2LnX2svJQHSPoo8A6Ky+oVwF/0Mlz6kCfpB8DewFbgXNs3t7ikpkm6GngdxeiZvwT+DrgBuAaYTDGi8Ntt99a5N6T0cS6/Af4JGA88Dqy0fUKramxGH+cxH9gN+HW52e22z25Jgduhj3N5MzANeJbiv6+zyz7MIau387D9rw3rH6Z4snOnj6absIiIiEq5DRUREZUSFhERUSlhERERlRIWERFRKWERERGVEhbRtsoROL/SsDxa0kZJ/1HDsZ4pRy/t/jlwAN8xTtIHd3ZtETtDrXNwR7TY74BDJI21vQU4HqjrOfot5eilO2Ic8EHg89uzk6RRtttmaJcYmnJlEe3uRuDE8vNc4OruFZKOlPQjSSsk/VDStLL9ryUtKj8fWs4R8oLtPbCkUeX8D8vKeRPeX7bvKelmSXdIWlUOBAfwSWBKeWWyQNLrGq+CJH2uHGoDSQ9L+pSkO4BTJU2R9B1JyyX9QNLLyu1OLeu/U9JwGh49hphcWUS7WwxcVP6jOwNYBBxTrrsXOMZ2l6TjgE9QjA/0j8B3JZ0MXAC8v4mhLcZKWll+fsj2ycCZFKMYv1LSbsBtkm4CHgFOtr2pnKjmdklLKObrOKT7CkXS6yqO+WvbR5Tb3kzxBvL9ko6iuDo5lmKAuRNsrx8uExXF0JSwiLZm+ydl/8FciquMRnsBXy7HCDIwptzn2fIv+J8Al9u+rYlD9XYb6o3ADElvazjeVIrxez4h6TUUQ01MZGDDr38NiisVisEvv14MbwYUQ3IA3AZcKekaikEyIwYkYREjwRKKOUteRzF2VrePAbfYPrkMlO82rJsKbKYYtn6gBHzI9tJtGosgGg+8wvbWcjyf3qby7WLbW8U9t/ld+XsX4PHe+kxsn11eaZwILJf0Ctu/7rldRJX0WcRIsAj4qO1VPdr34vcd3u/ubpS0F3AZ8Bpg74Yrg+21FPiApDHl975U0h7lcR8tg+L1wAHl9r8F/lvD/j8FpkvarbyF9IbeDlLO+PaQpFPL40jSy8vPU2z/2PZFwEaG/4RY0SIJi2h7ttfZvqyXVZcCl0hawbZX2Z8FFtq+j6Lf4ZOS9pXUIemL23HoLwJ3U8wvcBdweXmcrwIdklYB76IcGrv8i/+2skN6ge1HKEbdvav8vaKfY70TOFPSncBqiiG3ARaUneh3AT8Eap1NLdpXRp2NiIhKubKIiIhKCYuIiKiUsIiIiEoJi4iIqJSwiIiISgmLiIiolLCIiIhKCYuIiKj0/wEaAYXLJla4dwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "# Plot the score as a function of alpha\n",
    "ax.scatter(model.cv_results_['param_max_features'].data,\n",
    "          np.sqrt(-1 * model.cv_results_['mean_test_score']))\n",
    "ax.scatter([model.best_params_['max_features']], np.sqrt([-1*model.best_score_]), marker='o', color='r', s=40)\n",
    "ax.set_xlabel('Max. Features')\n",
    "ax.set_ylabel('RMSE (eV/atom)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save our best model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = model.best_estimator_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 3: Cross-validation Test\n",
    "Quantify the performance of this model using 10-fold cross-validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv_prediction = cross_val_predict(model, data[feature_labels], data['delta_e'], cv=KFold(10, shuffle=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute aggregate statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "r2_score 0.8236257047458013\n",
      "mean_absolute_error 0.190188108974116\n",
      "mean_squared_error 0.08647767302997238\n"
     ]
    }
   ],
   "source": [
    "for scorer in ['r2_score', 'mean_absolute_error', 'mean_squared_error']:\n",
    "    score = getattr(metrics,scorer)(data['delta_e'], cv_prediction)\n",
    "    print(scorer, score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,\n",
       "           max_features=12, max_leaf_nodes=None, min_impurity_decrease=0.0,\n",
       "           min_impurity_split=None, min_samples_leaf=1,\n",
       "           min_samples_split=2, min_weight_fraction_leaf=0.0,\n",
       "           n_estimators=20, n_jobs=-1, oob_score=False, random_state=None,\n",
       "           verbose=0, warm_start=False)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the individual predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO4AAADQCAYAAAAAnl3/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztnXl4VNX5+D9vQhYIS1gStrDLGkLYl7ogEBFUZFML2lorypcWK2i1UtuKys/Wte5W0WJVLKgUUQFRQRQXQEGiEUQWIZiwJKwBkhCSvL8/ZjK5MySTmWQms+R8nuc+3OXcc9+Z8M577jnvIqqKwWAILSICLYDBYPAeo7gGQwhiFNdgCEGM4hoMIYhRXIMhBDGKazCEIEZxDYYQxCiuwRCCGMU1GEKQeoEWoCa0aNFCO3bsGGgxDGFEdnY2Bw8eBDisqgmBlqcyQlpxO3bsyKZNmwIthiFMePvtt5k0aRIzZszg+eefzwy0PO4wQ2WDwc748eN55ZVXePbZZwMtSpUYxTXUeRYsWMC+ffuIiIjg+uuvJyIi+NUi+CU0GPzIc889x7Rp03j00UcDLYpXGMU11FkWLFjAzJkzGTdunFFcgyEUWLhwITfddBNjxozhrbfeIjo6OtAieYVRXEOdo6SkhKeeeoqLL76YpUuXEhMTE2iRvCakl4MMhuoQGRnJhx9+SFRUFPXr1w+0ONUiaCyuiLQTkbUisk1EtorIrEDLZAgv3n//fSZOnEhBQQHx8fHExcUFWqRqEzSKCxQDf1TVXsBQYKaI9AqwTIYwYc2aNUycOJHMzEzOnDkTaHFqTNAorqoeUNVv7PsngR+AtoGVyhAOrFu3jnHjxtGtWzc++ugj4uPjAy1SjQkaxbUiIh2BfsDGCq5NF5FNIrIpNze3tkUzhBgbNmzg8ssvp0OHDqxevZrmzZsHWiSfEHSKKyINgf8Bs1U1z/W6qs5X1YGqOjAhIWh9wA1BQkxMDL1792bNmjUkJiYGWhyfEVSzyiIShU1pX1fVpYGWxxC65OTkkJiYSL9+/fjyyy8RkXPaFBY7H8cGlTa4J2gsrti+2X8DP6jqPwMtjyF02bZtG7179+axxx4DqFBpQ52gUVzgfODXwEgRSbdvlwVaKENosXPnTkaNGkVkZCTjxo1z2za2nvMWSgSNuKr6ORB+P42GWmPPnj2MHDmS4uJiPv30U7p16xZokfxG0CiuwVATCgsLSUtL4/Tp06xdu5ZevcLbBcAoriEsiI2NZd68eXTr1o3U1NRAi+N3jOIaQpqcnBy2bt3KiBEjuPbaawMtTq1hFNcQshw5coS0tDSysrLYs2cPTZo0CbRItYZRXENIcvz4cUaPHs2OHTtYvnx5nVJaMIprCEHy8vIYM2YMGRkZvPPOO6SlpQVapFrHKK4h5HjppZfYvHkzS5YsYezYsfyUW+B0vXNCaMbYeoNRXEPIMXv2bIYPH86AAQMCLUrACCbPKYOPKSwu30KdE6fPcOO0m9n2426KSiOclLZN0/pO20+5BY7NHdbvJ9S+I6O4hqDn7Nmz/GrqNby84CU2rP8y0OIEBWaoHMaEmv+tK4XFUFxczG9+fR3L33uXOfc/Sv8R49menUde0VlHu8GdnGNsre+47t5/Q/n7MRbXELSUlJRw8403sHTJW9z+1wf45fU3B1qkoCGEf3MM4U5+fj67d+/ivnkPcOWvbwm0OEGFUVxDreNuImj/sQJUlaKiInomxfPZuk+JiYlh+daDTu3Surdy7H+154jTtYu6lg+d2zQNz6UhM1Q2BBWqygP33MW0aydy5syZkExWXhsYi2vwC+6samWTQqrK84/N4+UXnmHmH2aRdaIEEdvk0hXJrZzaerp8480EVCgtCRmLawga7r//fh59+EFunj6DRx57PCxTzvgKo7iGoODJJ5/k3nvv5dfX38ATTz9rlLYKvB4qi0gcUKiqJX6QxxAm7D9Wvn7qOkFkHZIu2rIPgJMtkxl1zQ1cdut9LN92AIBuTRs52qVnnnXqo0fbxo5962SUa//eDJVDaV23SosrIhEicq2IrBCRHGA7cMBe4+cRETnP/2IawpWftqajqrTt0o3r7riPiMjIQIsUEnjyG7MWWA38GfheVUsBRKQZMAJ4SETeVtWFNRFERBYAVwA5qtq7Jn0ZzsVdDuHqWih3uPNe2pR9FIC1yxbz/H13cOM9j3LBFVcD0CCqXHG35Bx37PdLdC4bUpn84GztwzVSyJM/U5qqnnU9qapHsSUv/589kXlN+Q/wDPCqD/oyBDmfv/82L9x/J72HDmfI6CsDLU7IUaXilimtiAwE/gJ0sN8ntsvapyLF9hZVXWevGWQIczZ8tJxn/zabXgOHccvD84mKNmu13uLNwOh14E4gAyj1jzhVIyLTgekA7du3D5QYIYe7IbCnw2NPh9vgPFz939b9jv2TR3N55q+z6JTcjxv//iJNGznXqM08XujYT0ksn5yyTka5Pi+UJpV8hTcfOVdV3/WbJB6iqvOB+QADBw7UAItj8JJGzRKY/uB8Oib3I6Z+g0CLE7J4o7hzReQlYA3gqAxsinPVHdxZtu3ZzoUVP9pTXgI1r7CUPenrKSo4TfdhafQcfFGl/Zyf1Myxn3Uqv9Jnu/Nysi4/hXJhL3d48zF+C/QAoigfKitgFNfgln3ff82b986geVInug4ege2/kKEmeKO4g1S1u78EEZFFwMVACxHJAuaq6r/99TxD7ZC5bQuL75lO44RWTJn3klmn9RHeKO6XItJLVbf5QxBVneqPfkMRb4Z3/p6kcde/9dqr6dlO146fLuLwnh/44O83Uy8unj43Pc6POQo5ufTo2NTRLiHO2fr27WCZhMqs+FlV4W7SzB2hFGTgzZ96KJAuInuwveM6loP8Ipkh5Nm36WOiGzRi4MyniI1PCLQ4YYU3ijvGb1IYnAjkBIqr1bFOOlnzPAGs21vu2ZR99LRjX1U5cOAUDQdMIbnnZdRvmuh036mC8n42/3zK6VqruPKgeGsuqep+J3XWV7kMVc2saPOncIbQ4+Shfaye9xvyc/chIkTF1a3SILWFV78xIpIKXGg//ExVv/W9SIZQ5VRuNp88/DtKzxahpSZ4zJ94rLgiMgu4mfLln4UiMl9Vn/aLZAaP8IXXkxXX/E1J8eVOEou+O+B0bcP35cencg+w9eXbKS4sIPmGR2neobwafFzDaKf7WsaXr7OmtHJ2wrA+zzpMd/WcCqVhrT/w5uNPA4ao6mkAEXkIWA8Yxa3jFJ08yrZX7qS4II9e1z9CXKsugRYp7PFGcQWwjn9K7OcMQYqnS0XulkGsfsYbXTIt5h6wheeVFBUQHd+GdmNnE9m0AwX5hTRuEuto166Fsz+y1coObNvM6ZqnWRnTM8utsdMSUh3BG8VdAGwUkbftxxMA4yBRhykpPIlE1CMyuj6dr5obaHHqFN4obitsbo8X2I9/C0zxuUQGn2G1slYLBc5Wyvpem5Fz0qndig1Zjv39mYcc+yVnTrN/+f8jMrYR7SfMJa6R87tqXIPy99pmDZzfca/p27ZSuawW110QfF20sla8SRZ3iap+o6pP2bctwFh/CWYIXkqLCti/4u8UHtlHs75XmMRuAaBKiysivwN+D3QWke8slxoBX/hLMENwUnq2kOyVD1KY+xNJl91Jo04DAy1SncSTofJ/gfeBfwBzLOdP2tPXGIIEd5XpGjeIqrTt7uPlXk/rdhx2anf4YLl3VHHecXK++DeFh34k4RfTSEguD89r3Mg5i0WiZcknxiUYyDo8dvXGquvLPJ7iSeqaE8AJwAQBGGiacgUNkvrQsMOgQItSp/HWc6op0BVwzPWr6jpfC2WoHHdW1RXXiZ/KWLqxPLJne8bPTtdKDmehpSWc2reZZkMmEx3fjLj2PQGQiPJ321YJzks+j08qT9Tp6tRhtf51fZKpunjjOXUTMAtIAtKxRQutB0b6RzRDMKClJeR+/Tqns7+lfrueNGibHGiRDHg3qzwLGARkquoIoB9w3P0thlBGS0vJ3byY09nf0izlSqO0QYQ3Q+VCVS0UEUQkRlW3i4jfMmLUZdz5FbvzLHK9lpdfPvFjTS4O8NWe8mH0z/vKh7Ilh/YAoFrK4fSlnP75G5oPmkKz/hOJT3BOSp5g8Yga0ME5Cujd78uH31f2bktlhGtOKH/jzdeUJSLxwDLgIxE5hlOOAkM4UXTiAKd+3kKz/pNp1n9ioMUxuODJOu4wYIOqlv317hWRtUATYJU/haureJMD2WqxrLmMwdnKrt3uPEGUmXXCsZ/3/abyC8U2Kx0T3Yi2Q28kru8lDgeLmBjnh08aUm5JXUuEWJd53EUcuU6uWSfffFE+JFwtuicf43rgWRHZgU1RV6nqp/4VyxAIVJVjuz4humEiDVv1IrphgvGKClI8Wcf9HYCI9MDm4vgfEWmCrRjYKuALU3IzPDi+ax0n9nxJo6T+NGzVK9DiGNzg8cBBVbdjK7H5uIjUx1ap72rgn4BP/N5EZAzwJBAJvKSqD/qi31CguiF41uFx1vF8p2urt5V7QWUfcA4eyNpjCdGLiOT47s84vusTGna9iOYXTgexLTg0bl4+6XTNRR2d+rAOj10D3SuTEdxPsPm6ul64DI1d8aQ+7rMicoH1nKoWqOpKVf2DqvpKaSOBZ7FZ9V7AVBExP/u1wIk96zm2/SPi2qTQ4oLpiHizSmgIBJ78Hu0AHhGR1sCbwCJ7ZJCvGQzsUtWfAERkMTAe8Ese52DDUyvrar225ZRPMlkLZgEUFZfXZjt10vmaHikPkC9BadCuHwkX3EyHbklO7Vq2LF/ysZYHgXP9n61YP48/atSaol9VoKpPAk+KSAds8bcL7EPlRdiUeIePZGkLWP3tsoAhro1MtT7fUVp8hoh6MTTtOxFUkQhjaUMFb9OzPqSq/bAFHEwAfvCbZJXLMV9VB6rqwIQEk2S7upz6+RuyPvgHZ0/mIiJGaUMMb3yV62F7/5wCjAI+Ae71oSzZQDvLcZL9XEjjiyGdNduh6/DUOjxesWm/07Xc3PJk46dPlO/n7VpP7qZFxLboTGRcM9r0OM9xrX9P5+TlKW3Kh8qeDo1d8cdaal0cHlvxxAHjEmwW9jLgK2AxML0s26MP+RroKiKdsCnsFOBaHz+jznNyz9fs/+gJYpp1oOX5NxNRL7rqmwxBhye/W3/GFkz/R1U95i9BVLVYRG4BPsC2HLRAVbf663m1RWWWwZ0Vcg3ds4a+PfLJbqdreYXlE1BR0c7D3fyT5ctDZ84UUXhoBwdWPUZM8450nfp3ImNsHkxN48szMt7yi47O/ec7B7pbcTfpVNcnj/yNJ5NTIwHExq+Azqp6v4i0B1qp6le+EkZVVwIrfdWfwZno5h1o3HMU8f0nO5TWEJp4MyPxHDCM8kwYJ7GtuxqCnIJDuygtyieiXgzNh/2GyJiGgRbJUEO8GcQMUdX+IrIFQFWPiYh5QfICd4nH3TnXr9tZ7qTfwqWe7A/7y99eDh509o5SLaXg4A72LbuP1qkX0e839zquJSU2cuxP7de6UrmsE1LerMea4bF/8cbinrV7NymAiCQApe5vMQSSgpyf2PfOPCIbxNNj/MxAi2PwId78Lj4FvA0kisgDwFXAX/0iVZjizgpZ/Xddc0VZszBuP+g8cRUfVz7oady4vI9TB3aT/d48ouMaM/gPz9CynbNH1P8NKXdesfoZu3pmVdfryZ23l/WzGstcPbwJMnhdRDZjW8MVYIKq1roDhqFqVJWtb/yDiKgYBv3+Seo3bRVokQw+xpN1XFFVBacIoUrbGJypbFnE3XJQ1innKJ+4qEjHfl6B8/LMzszyd9zO7cojeVrc+U9S2zSgVYfOAHSJd87CaMXq4OGrrIv+9lWu63jyjrtWRP5gX/5xICLRIjJSRF4BfuMf8QzecDI3m01vPoOWltI4McmhtIbww5Oh8hjgRmCR3avpOFAfm9J/CDzhp2ghgxcUHj/EipdmU3T6JN1HTKRRQuUJ2gyhjycOGIXY1nCfE5EooAVQoKomNasHVDb54nreuhyUf9Y5ocgGy2RVznHniZ7SUuVM3hG2zL+N0vwTXPePV2jT1eZ77C7Q3TphZJ0sCtccTeGGV38WVT0LHPCTLIZqUHTqGOkv3UZR3hF+/eB/aNMtJdAiGWoB83vqY9w5WbjD6hN82sXiHj9d5Nhv3ti5uFZSYj5bCo5zy2ML6NpvsNM1q5V1tZyhsCRj/J0rx+uvQ0S6+TB43lBNSktLiIiIpFv/Ydy35DPqxzWq+iZD2FCd6Omb7TVzDQGiqOA0y++/kR2fvgtglLYOUp0ByFFghj1d67dAuqp+41uxQpfq5o6yDmtdy4WM6dXCsb9572EWPXArOTszuPDqm0ix+BwP7tTcY1n8PfT0RWJzMzyuHE8cMLoDOyxOGP8QkTXAj0Bf4ALAKG4tUFRYwJv3/Y6sbd8w4U+P0eP80YEWyRAgPPlNWwJ0sFcy+A7IsP8bY69oYKoa+ABrkSyrpxRAg6hISoqLeeRPM8j8dgOPPP0iE6+xJQexTmr5ys/YFxhvKf/iyTpuiojEAH2wBbmfBq4AkkUEVTWOsLVARGQk3VP7M3HCRIfSGuouHr1FqOoZ4GsROaWqfyg7b69Qb/AjJcXFHDjwM63bd2TKjNvpldik6psMYY+3r/9OgQT+zEEVCKpbBsTTSRSrM7+rJ9NAypONlw1/S0pKuOePM/h87Ye8vWYTzVokVDAErdzryd2asr8nfswarH/xZHLqWWyTT99gC+cz1AKlpaX8v7tnsXLZm9xy599o1sLkkDaU48lv4bfYZo+vBxqJyDZgK7bSINtU9Q0/ylereGpl3bVzDYK3WlZryJxrOyuN6tfj3jm3seyN17jl9jncduefK5TDVZZAWjbj41y7VOmAYa8c8AdVHa6qLYDRwAIgH9skVY0RkatFZKuIlIqIT4qIhTJvLHyZhS/P5+aZtzH7rr8FWhxDEOL176KqZgFZInIK8FXxnu+BScALPuovpJl49bWICNdcd4MpLG2oEK8UV0T6YasucA1wEOgB1DgLWVkKnNr+T+rN5I2n2Svc1Ym1ehNVVM7jjYUvc+nl4+nZrS09b/t95cIFIb4aGptJLc/wpD5uNxGZKyI/Ai8Ch4GLVXUINvfHWkVEpovIJhHZlJubW9uP9xsvPP1P7r59Jq++9K9Ai2IIATz5TduOra7PVaqa4XLN4zxTIrIaqMhZ4y+q+o6n/ajqfGA+wMCBA2uU58rdL7prGRArrksynloJayjd6h/LK8K/+9qLvPjQPVz9yyk8MO++ak/0hMMEUSjKHAg8+ZomYSvA9aFd+d4EVtmD6j1GVdOqIV/Y8/4br/DiQ/cwfuIk/v3yq0RGRlZ9k6HO44nL4zJgmYjEYasQPx14SURWAr5JCRgkWC2W1Tq6a+cNVl/iK5JbUVBQwOz/vsjll1/Om4sXER1deRnLyp5d3eUg8y4Z2niTV/k0tqp9/7W7Ol4NdPCFECIyEXgaSABWiEi6ql7qi76Dmfr16/PZZ5/RtGlToqNNNReD51SrDLmqHrOv7470hRCq+raqJqlqjKq2DHelXfXeMm6++WZKSkpo3bo1sbGxVd9kMFgI+0GSNxM27pZ83PVhbevueZ0T6vPee+8x6/+uZ+CgwRw/VUhcXFyNnl1dzPA4tKmWxTVUjw8++ICrrrqKfv36sey9lQ6lNRi8pUa/uyIyW1Wf8JUw/sBXEzaeWtXKrOXHH3/MhAkT6NWrFx988AH1G1UenuePyCRDeFFTi3u7T6SoA6gqffr04aOPPqJpUxPGbKgZNf29No60VXD06FGaNWvGqFGj2LBhg/E9NviEmlrckK7QV1jsvMXWK9/cYW0XW8/mZVW2Wc9v++4bunbtysKFCwFnX+z9xwqcNldZrJu7Z1fWzhDeeBJIfxKbglZkKhr4XKIwISMjg9GjR9OwYUMuvPDCWn9+x44dyczMrPXnhhEDRMRvhik2NvZQQUFBtfO1eeI5FXLZtj2dZHJ3nyvuashafZcLi+HH7dsZnZZGdEwMH3/8MR06dDinf39nQczMzMSULA5eRKRlTe73xOK+6+66ql5ZEwHCjWPHjjH20lGICKs+/JguXboEWiRDGOLJ5NQwIAubu+NGzISUW5o2bcodd85h+MUj6Na9e6DFMYQpUtVwSkQigUuAqdhyK68AFqnqVv+L554BAwbqFxs3AZWvpVZ3TdTba1lZWRw6dIgBAwZU+gxrH65hg74eOttzXvu0T4PvsP99qm0EPck5VaKqq1T1N8BQYBfwiYjcUt2HhhsHDhxg5MiRTJw4kTNnzgRanEoREX71q185jouLi0lISOCKK5xTh02YMIGhQ4c6nbv33ntp27Ytffv2dWzHj7uvbb5582ZSUlI477zzuPXWWyv8Idm+fTvDhg0jJiaGRx991Onak08+Se/evUlOTuaJJzz389m7dy9JSUmUlpY6ne/bty8bN24EbH+z0aMrL+Fy/PhxnnvuOY+fWeuoapUbEIMtLvctbEH1fwPaenKvP7f+/QdowVnVgrOqu3PyHZuVsutlm7trlbVzR05Ojvbs1Uvj4uL040+/8PxGL6iOnLY/rTNxcXGampqq+fm272jlypWampqql19+uaPNsWPHNCkpSXv06KG7d+92nJ87d64+8sgjXsk9aNAgXb9+vZaWluqYMWN05cqV57Q5dOiQfvXVV3r33Xc79Z+RkaHJycl6+vRpPXv2rI4aNUp37tzp8bOHDRumn3zyieP4hx9+0M6dOzuOFyxYoI8++mil9+/Zs0eTk5M9fp632P8+1f6/70nqmleB9UB/4D5VHaSq81Q1u4pbw56jR4+SlpbG3j17ePvdFQz7xS8CLVKVXHbZZaxYsQKARYsWMXXqVKfrS5cuZdy4cUyZMoXFixdX+zkHDhwgLy+PoUOHIiJcf/31LFu27Jx2iYmJDBo0iKgo51jkH374gSFDhtCgQQPq1avH8OHDWbp06Tn35+bmMnnyZAYNGsSgQYP44osvAJg6daqT/IsXL2bKlCmO41WrVjF27FhOnTrFqFGj6N+/PykpKbzzji0Zy5w5c9i9ezd9+/blzjvvRFW588476d27NykpKbzxhi0r8SeffMLw4cMZP348nTt3Zs6cObz++usMHjyYlJQUdu/eXe3v0C1VaTZQCpy0b3mW7SSQV5NfjZpugba4d911l8bExOjy9z/02lJ7gy8t7rfffquTJ0/WgoICTU1N1bVr1zpZ3LS0NF23bp3++OOP2rt3b8f5uXPnaps2bTQ1NVVTU1P14osvVlXV7OxsHTt27DnP+vrrr3XUqFGO43Xr1jk9xxVXi75t2zbt2rWrHj58WE+fPq1Dhw7VW2655Zz7pk6dqp999pmqqmZmZmqPHj1UVfXgwYPaqlUrPXvW9iX16NFDMzIyVFW1uLhYU1NTVVX17NmzeuLECVVVzc3N1S5dumhpaek5FnfJkiWalpamxcXFevDgQW3Xrp3u379f165dq02aNNH9+/drYWGhtmnTRu+55x5VVX3iiSd01qxZFX5eamhxPVnHDdoIIpHyCZ/KJne8Wce1ZqhwzYBR0STXvHnzmDRpEoMHD/ZI3urWjPVlIEGfPn3Yu3cvixYt4rLLLnO6dujQIXbu3MkFF1yAiBAVFcX3339P7969Abjtttu44447nO5p06YNK1eu9J2Adnr27Mldd93F6NGjiYuLo2/fvhWm9Vm9ejXbtm1zHOfl5XHq1ClatmxJ7969WbNmDS1btqRevXqOz7Fx40aGDBkC2AzX3Xffzbp164iIiCA7O5tDhw6d85zPP/+cqVOnEhkZScuWLRk+fDhff/01jRs3ZtCgQbRu3RqALl26ON6dU1JSWLt2rc+/GzBhfV6Tn5/PjBkzyMnJISoqymOlDSauvPJK7rjjjnOGyW+++SbHjh2jU6dOdOzY0aHg1aFt27ZkZWU5jrOysmjbtq1XfUybNo3Nmzezbt06mjZtSrdu3c5pU1payoYNG0hPTyc9PZ3s7GwaNmwIlA+XFy9e7PRZ33//fcaMGQPA66+/Tm5uLps3byY9PZ2WLVtSWFjolZwxMTGO/YiICMdxREQExcX+8UWtc4pr9fPdnp3ntFXWrsziFRYWMn78eF588UWWffCpwz/ZUzon1Hdsrnjqc+zOj9lTbrzxRubOnUtKSorT+UWLFrFq1Sr27t3L3r172bx5c7Xfc1u3bk3jxo3ZsGEDqsqrr77K+PHjveojJycHgH379rF06VKuvfbc8qKjR4/m6aefdhynp6c79idNmsTKlSt54403nN5v16xZQ1qaLXfhiRMnSExMJCoqirVr1zrcRBs1asTJkycd91x44YW88cYblJSUkJuby7p16wL6o22iOT2kqKiIyZMns2bNGh5+aj5pY3xSfSUgJCUlceuttzqd27t3L5mZmU7LQJ06daJJkyaOJZTHH3/cETABsGzZMqKjo7npppsqHC4/99xz3HDDDRQUFDB27FjGjh0LwPPPPw/AjBkzOHjwIAMHDiQvL4+IiAieeOIJtm3bRuPGjZk8eTJHjhwhKiqKZ599lvj4+HOe8dRTTzFz5kz69OlDcXExF110kaP/+Ph4hg0bxsGDB+ncuTNgm8yKjY2lUSObJ+91113HuHHjSElJYeDAgfTo0QOA5s2bc/7559O7d2/Gjh3Lww8/zPr160lNTUVEePjhh2nVqhXbt2+v3h+hhlTpgBHMDBw4UDdt2lTt+10Lb1mrC1it4tmzZ7nmmmtYtmwZ8+fPZ9SEXznd5wvnCV+X+DQOGBWzcOFCsrKymDNnTkDlqKkDRlBYXBF5BBgHFAG7gd+qqvvVfeBMcaljqGqdTPI00blrGZDKFPDYsWN88+33zPrbg3S/eJKjfm053qdydVfixB3eKLXhXKwOKKFMUCgu8BHwZ1UtFpGHgD8DdwVYJkpLS1FVEhMTeemdtcTE+jeix2DwlKBQXFX90HK4AbjKk/ti6kVUOUx1tUKeDmt3HTrNX+64hdOnTvH4v17mwp5JHt3nSm3mhDL5p+oOwTirfCPwfiAFUFXuv/uPvLnwP3Ts3MX+n+GVAAAMKklEQVSUBTEEHbX2G+1J0S8R+QtQDLzupp/p2Mqg0K59e4dF9dTaLN9aXmwrrfu54qgqc/50B68teIGbfj+b2+6655w2/rBspiSIwRtq7b+IVlH0S0RuwFbhfpS6mQ5VS7W+AQNqVq2vIv7xwDyeeuKf/HraDObMfaDOJHdbtmwZK1asIC8vj2nTprmNnDEEnqAYKovIGOBPwJWqmh9IWYZfPJJbbp3NPQ88GpZK+8ILL9CqVStSU1Pp0qULr776KmAL5XvxxRd5/vnnHQ703rBq1Sq6d+/Oeeedx4MPPlhhm8cff5zk5GR69+7N1KlTHR5KP//8MyNGjKBXr14kJyfz5JNPVv8D1hVq4ujsqw1bjO/PQLp9e96T+6xBBp4GCFTWbsuWLZW28yZ4wNOgBn9DBUEGqqozZ87Uf/3rX6qqunHjRm3evLnT9dtvv103b97s1bOKi4u1c+fOunv3bj1z5oz26dNHt27d6tQmKytLO3bs6AgpvPrqq/Xll19WVdX9+/c7npmXl6ddu3Y95/5wA3+H9dUGqnqeqrZT1b72bUZtPv+pp56iX79+jnC3cOa7776juz2lTqdOnRxVAlWVu+66i7Fjx9K/f3+v+vzqq68477zz6Ny5M9HR0UyZMsURHmeluLiYgoICiouLyc/Pp02bNoDNPbLsmY0aNaJnz55kZ9f5qFG3BIXiVpey6KCK/Iorw7XdCy+8wKxZs5g4caLTe503fVbWf0VOFt7252syMjLo3r07qsozzzzDAw88AMDTTz/N6tWrWbJkicNlEGzxu/v373fbZ3Z2Nu3atXMcJyUlnaN4bdu25Y477qB9+/a0bt2aJk2aVPgevXfvXrZs2eKI3jFUQk3MdaC3AQMG1GS0ogsWLFBAL7/8cj1z5kyN+go2qGCovG/fPo2MjNTU1FRt0aKFjhw5UktLS2v8rLfeekunTZvmOH711Vd15syZTm2OHj2qI0aM0JycHC0qKtLx48fra6+95tTm5MmT2r9/f/3f//5XY5mCHcJhqBwIfvzxR2666SYuueQSlixZUicKS2dkZHDRRReRnp7Ojh072L59O+vXr69xv23btuXnn392HFcUwrd69Wo6depEQkICUVFRTJo0iS+//NJx/ezZs0yePJnrrruOSZMm1VimcCekFVe16pId6Zl5TlsZ3bt3Z+GiN1m8ZBnUi/XKz9cXoXWB4LvvvqNfv36ALY3stdde65P3+kGDBrFz50727NlDUVERixcv5sorndNtt2/fng0bNpCfn4+qsmbNGnr27AnYRn3Tpk2jZ8+e3H67qSPnCSGtuNVh+fLljrxEEydNpkGDulNFJSMjw6G4AOPGjasye4Un77j16tXjmWee4dJLL6Vnz55cc801JCcnO90/ZMgQrrrqKkdup9LSUqZPnw7AF198wWuvvcbHH3/syCDpj6wa4URIh/VZ8ypbqSwH8kcffsBVE6+k36ChvL70fbokOittdb2XfOH15GnZFE/7N2F9wY3f8yqHC59+spZrJk+gS9fuPLfgv2HpXGGoO9QJxf3i88+ZNP4KOnfpwitvLSe+abNAi2Qw1IiQdme3Znm04jph9NorL9M2KYkVq1bTMalVpe2C1bk/WOUyBI6w/i+hqogIz/zrBY4ePUpiYmKgRTIYfEJYDpVj68HOHzIYPfIijubsp2FsPdq3SXTryeTtZFR1KtlXJXNlfYTi0pPBv4Slxd2+fTtpaWnUq1ePggLP06eGEx06dDATcEFMbGzsuVnXvSDsFHfXrl2MHDkSEeHjj+tuYem9e/cGWoSQRkQ2q+rAQMtRGSGtuGWeU2VkZmZyyahRFBUV8cHqT+jQxT+FpWt7sshMThlcCat33NjYWDp06Mjy9z8i2V4nxmAIR0L6t7yopJT9xwo4euQwqV1a0aFtS9Z9+onf3+1CZRnJEL6EvMU9cjiXaydcym9/+1sAMyFjqBOEtOKWFJfwm6uvYF/mHqZNmxZocQyGWiOkB3mZe3ZRXFzMe++9x4gRI3zat7vhsBkaGwJNSP8XzM/P57333jOpRA11jpAO6xORXCDTj49oARz2Y//Bgvmc59JBVRP8KUxNCGnF9TcisimYF+F9hfmcoUdIT04ZDHUVo7gGQwhiFNc98wMtQC1hPmeIYd5xDYYQxFhcgyEEMYprMIQgRnGrQEQeEZHtIvKdiLwtIvGBlslXiMgYEflRRHaJyJxAy+MPRKSdiKwVkW0islVEZgVaJl9g3nGrQERGAx+rarGIPASgqncFWKwaIyKRwA7gEiAL+BqYqqrbAiqYjxGR1kBrVf1GRBoBm4EJof45jcWtAlX9UFXLPJc3AEmBlMeHDAZ2qepPqloELAbGB1gmn6OqB1T1G/v+SeAHoK37u4Ifo7jecSPwfqCF8BFtsRUTLyOLMPgP7Q4R6Qj0AzYGVpKaE9JBBr5CRFYDrSq49BdVfcfe5i9AMfB6bcpm8A0i0hD4HzBbVfOqah/sGMUFVDXN3XURuQG4Ahil4TMpkA20sxwn2c+FHSIShU1pX1fVpYGWxxeYyakqEJExwD+B4aqaG2h5fIWI1MM2OTUKm8J+DVyrqlsDKpiPEVtKlFeAo6o6O9Dy+AqjuFUgIruAGOCI/dQGVZ0RQJF8hohcBjwBRAILVPWBAIvkc0TkAuAzIAMotZ++W1VDuo6nUVyDIQQxs8oGQwhiFNdgCEGM4hoMIYhRXIMhBDGKazCEIEZxDYYQxCiuwRCCGMWtBiJSIiLp9vjOb0XkjyIS4XKtbBtg2T8oItmW4+gK+p4gIioiPVzO/5+I/Mvl3Pci0rMKWeuLyKf2ML7K2qwVkUtdzs0ue56IPC8i57u5P15Efu9ODm8QkWgRWWf37jJUgFHc6lGgqn1VNRlbPOtYYK7LtbJtc9k+8DzwuOVaUQV9TwU22f+1kgJ8U3YgIrFAR2xui+64EViqqiVu2iwCpricm2I/DzAUW0hjZcQDPlNc+/eyBvilr/oMN4zi1hBVzQGmA7dIDUsF2iNYLgZu4lzF7YNFcbEp8o4qFBLgOqAswulXIvKV3dq/YLHCS4DLy0YA9vC3NsBndovueI6ILBORzfbRxnT7/Q8CXez9PmJvd7t9RPC9iMwu69eeTeQ/IrJDRF4XkTQR+UJEdorIYIvcy+yyGypCVc3m5QacquDccaAlUAKk27e3XdrcC9zhpt/rsEWwgE1JB1iuHcFWbmWvfTsM/KcKOaOBg/b9nsB7QJT9+Dngekvb5cB4+/4c4FH7/u3AjZZ2zez/1ge+B5pjs/zfW9oMwOYbHAc0BLZii4PtiC00MgWb0dgMLAAEWxD/MksfkUBuoP/WwbqZdwjfU6C2YXF1mAq8aN9/0368WUTaYftP7HjvFZFngD2W4z9jU6KXVHW7/XQLbD8oYIsCGgB8bR8Y1AdyLM8uGy6/Y/+3rG7ppcBvLe1uFZGJ9v12QFfgoMvnuADbj9Zpu2xLgQuBd4E9qpphP78VWKOqKiIZ2BQbAFUtEZEiEWmktswVBgtmqOwDRKQzNkubU1VbN300A4YAq+yn3gR+aR9+p2CzWlZ6Ad/Z7x2CTcn3WpQWoACILXsE8IqWv193V9V7LW3fAUaJSH+ggapuFpEGQLyq7rc/52IgDRimqqnAFkv/nnLGsl9qOS7l3PjwGKDQy/7rBEZxa4iIJGCbdHpG7WO8anIVsFJVzwCo6k/AAWyWqg/gmtwsGdtwFOBH4FNVfcbaQFWPAZH2iaw1wFUikmiXu5mIdLC0PQWsxTZ0LZuUGmE/V0YT4Jiq5ttnvYfaz58EGlnafQZMEJEGIhIHTLSf8xgRaQ4cVtWz3txXVzBD5epRX0TSgShs72yvYQu2rwlTgVQR2Ws519x+vjG2d1DAYZ1FVcuGqH2Bbyvp90PgAlVdLSJ/BT60L12dBWbiXKZ0EfA25TPMY7FNXJWxCpghIj9g+7HYAKCqR+wTTN8D76vqnSLyH+Ar+30vqeoW+6SXp4wAVnjRvk5h4nHDAPus7eequqmCa/2B21T119Xo9xtgSCCsnv29eI6qVrXcVScxilsHEJEbsb3fVrV0FBTYl6WmqOqrgZYlWDGKazCEIGZyymAIQYziGgwhiFFcgyEEMYprMIQgRnENhhDEKK7BEIIYxTUYQpD/D/EMnOKy4r9nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 216x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.hist2d(pd.to_numeric(data['delta_e']), cv_prediction, norm=LogNorm(), bins=64, cmap='Blues', alpha=0.9)\n",
    "\n",
    "ax.set_xlim(ax.get_ylim())\n",
    "ax.set_ylim(ax.get_xlim())\n",
    "\n",
    "mae = metrics.mean_absolute_error(data['delta_e'], cv_prediction)\n",
    "r2 = metrics.r2_score(data['delta_e'], cv_prediction)\n",
    "ax.text(0.5, 0.1, 'MAE: {:.2f} eV/atom\\n$R^2$:  {:.2f}'.format(mae, r2),\n",
    "        transform=ax.transAxes,\n",
    "       bbox={'facecolor': 'w', 'edgecolor': 'k'})\n",
    "\n",
    "ax.plot(ax.get_xlim(), ax.get_xlim(), 'k--')\n",
    "\n",
    "ax.set_xlabel('DFT $\\Delta H_f$ (eV/atom)')\n",
    "ax.set_ylabel('ML $\\Delta H_f$ (eV/atom)')\n",
    "\n",
    "fig.set_size_inches(3, 3)\n",
    "fig.tight_layout()\n",
    "fig.savefig('oqmd_cv.png', dpi=320)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
